# vrptw_deadline_split_ga_tabu_energy_capacity.py
import math, random, copy, json
from dataclasses import dataclass, replace
from typing import List
import numpy as np
import pandas as pd

# ====== Numba: 优雅降级 ======
try:
    from numba import njit
except Exception:
    def njit(*args, **kwargs):
        def wrap(f): return f
        return wrap

# =======================
# 数据结构
# =======================
@dataclass
class Customer:
    id: int
    x: float
    y: float
    service: float         # 服务时间（分钟）
    deadline: float        # R_i: 本单“自出仓起”必须在多少分钟内送达
    weight: float = 1.0    # 迟到单价/权重（元/分钟）
    demand: float = 0.0    # 该单重量（kg）

@dataclass
class Depot:
    x: float
    y: float

@dataclass
class Instance:
    customers: List[Customer]
    depot: Depot
    speed: float = 40.0/60.0   # km/min
    work_budget: float = 3200.0

@dataclass
class EvalParams:
    # 目标权重
    late_weight: float = 1.0
    lambda_duration: float = 0.02
    overtime_penalty: float = 1e9
    seed: int = 123
    # ===== 资源与能耗参数 =====
    max_total_load_kg: float = 500.0   # 最大载重（含自重）
    curb_weight_kg: float = 100.0      # 车辆自重
    max_payload_kg: float = 400.0      # 可装货物上限（= 500-100）
    energy_max_range_km: float = 100.0 # 最大载重下的最远续航（km）

@dataclass
class Solution:
    order: List[int]  # customers 的 id 排列（1..N）

@dataclass
class DecodeResult:
    cost: float
    total_late_cost: float
    total_duration: float
    total_energy: float
    routes: List[List[int]]
    route_energies: List[float] = None  # 每趟的电量

# =======================
# 工具
# =======================
def seed_all(s=123):
    random.seed(s)
    np.random.seed(s)

def euclid(ax, ay, bx, by) -> float:
    return math.hypot(ax - bx, ay - by)

class TravelTime:
    """ 预计算（时间矩阵 min & 距离矩阵 km）"""
    def __init__(self, ins: Instance):
        n = len(ins.customers)
        self.n = n
        d = ins.depot
        sp = ins.speed  # km/min

        xs = np.array([c.x for c in ins.customers], dtype=np.float64)
        ys = np.array([c.y for c in ins.customers], dtype=np.float64)

        self.service   = np.array([c.service  for c in ins.customers], dtype=np.float64)
        self.deadline  = np.array([c.deadline for c in ins.customers], dtype=np.float64)
        self.penalty_w = np.array([c.weight   for c in ins.customers], dtype=np.float64)
        self.demand    = np.array([c.demand   for c in ins.customers], dtype=np.float64)

        # depot -> customer 距离（km）与时间（min）
        dx = xs - d.x; dy = ys - d.y
        dep2c_dist = np.sqrt(dx*dx + dy*dy)          # km
        self.dep2c      = dep2c_dist / sp            # min
        self.dep2c_dist = dep2c_dist                 # km
        self.c2dep      = self.dep2c.copy()          # min
        self.c2dep_dist = self.dep2c_dist.copy()     # km

        # customer -> customer
        self.c2c      = np.zeros((n, n), dtype=np.float64)  # min
        self.c2c_dist = np.zeros((n, n), dtype=np.float64)  # km
        for i in range(n):
            dx_i = xs - xs[i]; dy_i = ys - ys[i]
            dist_row = np.sqrt(dx_i*dx_i + dy_i*dy_i)   # km
            self.c2c_dist[i, :] = dist_row
            self.c2c[i, :]      = dist_row / sp        # min

# =======================
# Numba JIT: 区间表 + 一维 DP（含容量+能量）
# =======================
@njit
def _build_seg_tables(order_idx,
                      dep2c, c2dep, c2c, service, deadline, penalty_w,
                      demand, dep2c_dist, c2c_dist, c2dep_dist,
                      late_w, lambda_dur,
                      max_payload_kg, max_total_load_kg, curb_weight_kg, energy_max_range_km):
    n = len(order_idx)
    INF = 1e18
    seg_cost = np.full((n, n), INF, dtype=np.float64)
    seg_dur  = np.full((n, n), INF, dtype=np.float64)
    seg_energy = np.zeros((n, n), dtype=np.float64)

    # 电池100%容量
    denom_energy = max_total_load_kg * energy_max_range_km

    for i in range(n):
        idx_i = order_idx[i]

        # ---- 时间/迟到初始化 ----
        cumul = dep2c[idx_i]
        late_sum = max(0.0, cumul - deadline[idx_i]) * penalty_w[idx_i]
        cumul += service[idx_i]
        prev = idx_i

        # ---- 载重 & 能耗簿记 ----
        # 路径距离（km）：depot -> ... -> last
        path_km_to_last = dep2c_dist[idx_i]
        # 起点载重（kg）：全部订单重量
        cap_sum = demand[idx_i]
        # 货物能耗项： sum(demand[j] * distance_from_depot_to_j)
        energy_payload_num = demand[idx_i] * path_km_to_last

        for j in range(i, n):
            if j > i:
                idx_j = order_idx[j]
                # 时间/迟到
                travel_t = c2c[prev, idx_j]
                arrive   = cumul + travel_t
                late_sum += max(0.0, arrive - deadline[idx_j]) * penalty_w[idx_j]
                cumul    = arrive + service[idx_j]

                # 距离/载重/能耗更新
                path_km_to_last += c2c_dist[prev, idx_j]         # km
                cap_sum += demand[idx_j]                          # kg
                energy_payload_num += demand[idx_j] * path_km_to_last
                prev = idx_j

            # 收车回仓
            back_t   = c2dep[prev]
            back_km  = c2dep_dist[prev]
            dur = cumul + back_t
            seg_dur[i, j] = dur

            # 能耗计算 - 检查是否超出100%电量
            route_km = path_km_to_last + back_km
            energy_total_num = curb_weight_kg * route_km + energy_payload_num
            energy_frac = energy_total_num / denom_energy

            # 可行性检查：容量 + 电量都不能超限
            if (cap_sum <= max_payload_kg) and (energy_frac <= 1.0):
                seg_cost[i, j] = late_w * late_sum + lambda_dur * dur
                seg_energy[i, j] = energy_frac
            else:
                seg_cost[i, j] = INF
                # 继续延伸只会更重/更耗电，提前停止该 i 的扩展
                break

    return seg_cost, seg_dur, seg_energy

@njit
def _dp_from_seg(seg_cost):
    n = seg_cost.shape[0]
    dp   = np.empty(n, dtype=np.float64)
    prev = np.empty(n, dtype=np.int64)
    for j in range(n):
        dp[j] = seg_cost[0, j]
        prev[j] = -1
        for i in range(1, j+1):
            cand = dp[i-1] + seg_cost[i, j]
            if cand < dp[j]:
                dp[j] = cand
                prev[j] = i-1
    return dp, prev

# =======================
# DP 精评（调用 JIT）
# =======================
def split_dp(order_ids: List[int], id2idx_map: dict, tt: TravelTime, ep: EvalParams, work_budget: float) -> DecodeResult:
    n = len(order_ids)
    order_idx = np.array([id2idx_map[cid] for cid in order_ids], dtype=np.int64)

    seg_cost, seg_dur, seg_energy = _build_seg_tables(
        order_idx,
        tt.dep2c, tt.c2dep, tt.c2c, tt.service, tt.deadline, tt.penalty_w,
        tt.demand, tt.dep2c_dist, tt.c2c_dist, tt.c2dep_dist,
        ep.late_weight, ep.lambda_duration,
        ep.max_payload_kg, ep.max_total_load_kg, ep.curb_weight_kg, ep.energy_max_range_km
    )
    dp, prev = _dp_from_seg(seg_cost)

    # 回溯
    routes = []
    route_energies = []
    j = n - 1
    while j >= 0:
        i = (prev[j] + 1) if prev[j] != -1 else 0
        seg = order_ids[i:j+1]
        routes.append(seg)
        route_energies.append(seg_energy[i, j])
        j = prev[j]
    routes.reverse()
    route_energies.reverse()

    # 统计
    total_late_cost = 0.0
    total_duration  = 0.0
    total_energy   = 0.0
    s = 0
    for r in routes:
        i = s
        j = s + len(r) - 1
        late_only = seg_cost[i, j] - ep.lambda_duration * seg_dur[i, j]
        total_late_cost += late_only
        total_duration  += seg_dur[i, j]
        total_duration += seg_energy[i, j]/0.04  #每分钟充电0.25%
        total_energy   += seg_energy[i, j]
        s = j + 1

    overtime = max(0.0, total_duration - work_budget)
    total_cost = dp[-1] + ep.overtime_penalty * overtime
    return DecodeResult(total_cost, total_late_cost, total_duration, total_energy, routes, route_energies)

# =======================
# 贪心快评（含容量+能量）
# =======================
def split_greedy(order_ids: List[int], id2idx_map: dict, tt: TravelTime, ep: EvalParams, work_budget: float) -> DecodeResult:
    INF = 1e18
    routes = []
    route_energies = []
    total_late_cost = 0.0
    total_duration = 0.0
    total_energy = 0.0

    n = len(order_ids)
    if n == 0:
        return DecodeResult(0.0, 0.0, 0.0, 0.0, [], [])

    # 电池100%容量
    denom_energy = ep.max_total_load_kg * ep.energy_max_range_km
    demand = tt.demand

    def eval_single(idx):
        # 单点一趟：容量 & 能量可行性
        payload = demand[idx]
        if payload > ep.max_payload_kg:
            return INF, tt.dep2c[idx] + tt.c2dep[idx], 0.0
        path_km = tt.dep2c_dist[idx] + tt.c2dep_dist[idx]
        energy_num = ep.curb_weight_kg * path_km + payload * tt.dep2c_dist[idx]
        energy_frac = energy_num / denom_energy
        # 检查是否超出100%电量
        if energy_frac > 1.0:
            return INF, tt.dep2c[idx] + tt.c2dep[idx], 0.0
        # 成本
        cumul = tt.dep2c[idx]
        late = max(0.0, cumul - tt.deadline[idx]) * tt.penalty_w[idx]
        cumul += tt.service[idx]
        # 添加充电时间（每趟用电量/0.04）
        charge_time = energy_frac / 0.04
        dur = cumul + tt.c2dep[idx] + charge_time
        return ep.late_weight * late + ep.lambda_duration * dur, dur, energy_frac

    cur_route = []
    cur_late_sum = 0.0
    cur_cumul = 0.0
    prev = -1
    cur_duration_if_close = 0.0

    # 载重/能耗簿记（当前趟）
    cur_payload_kg = 0.0
    cur_path_km_to_last = 0.0
    cur_energy_payload_num = 0.0

    for cid in order_ids:
        idx = id2idx_map[cid]
        if not cur_route:
            # 开新趟
            cur_route = [cid]
            cur_payload_kg = demand[idx]
            cur_path_km_to_last = tt.dep2c_dist[idx]
            cur_energy_payload_num = demand[idx] * cur_path_km_to_last

            cur_cumul = tt.dep2c[idx]
            cur_late_sum = max(0.0, cur_cumul - tt.deadline[idx]) * tt.penalty_w[idx]
            cur_cumul += tt.service[idx]
            prev = idx
            
            # 计算当前趟的电量和充电时间
            route_km = cur_path_km_to_last + tt.c2dep_dist[prev]
            energy_num = ep.curb_weight_kg * route_km + cur_energy_payload_num
            energy_frac = energy_num / denom_energy
            charge_time = energy_frac / 0.04  # 充电时间
            cur_duration_if_close = cur_cumul + tt.c2dep[prev] + charge_time

            # 可行性：容量+能量（100%电量限制）
            if (cur_payload_kg > ep.max_payload_kg) or (energy_frac > 1.0):
                return DecodeResult(INF, INF, 1e9, INF, [])
            continue

        # 方案 A：并入当前趟
        travel_t = tt.c2c[prev, idx]
        arrive = cur_cumul + travel_t
        add_late = max(0.0, arrive - tt.deadline[idx]) * tt.penalty_w[idx]
        cumul_new = arrive + tt.service[idx]
        
        # 计算并入后的电量和充电时间
        new_payload = cur_payload_kg + demand[idx]
        new_path_km_to_last = cur_path_km_to_last + tt.c2c_dist[prev, idx]
        new_energy_payload_num = cur_energy_payload_num + demand[idx] * new_path_km_to_last
        route_km_with_back = new_path_km_to_last + tt.c2dep_dist[idx]
        energy_num_new = ep.curb_weight_kg * route_km_with_back + new_energy_payload_num
        energy_frac_new = energy_num_new / denom_energy
        charge_time_A = energy_frac_new / 0.04  # 充电时间
        
        dur_if_close_A = cumul_new + tt.c2dep[idx] + charge_time_A
        delta_A_late = add_late
        delta_A_dur  = dur_if_close_A - cur_duration_if_close
        cost_A = ep.late_weight * delta_A_late + ep.lambda_duration * delta_A_dur

        violate_capacity = new_payload > ep.max_payload_kg
        violate_energy = energy_frac_new > 1.0  # 100%电量限制
        violate_A_budget = (total_duration + (dur_if_close_A - cur_duration_if_close)) > work_budget + 1e-9
        violate_A = violate_capacity or violate_energy or violate_A_budget

        # 方案 B：单开一趟
        cost_B_single, dur_B_single, energy_B_single = eval_single(idx)

        if (cost_A <= cost_B_single) and (not violate_A):
            # 接受 A
            cur_route.append(cid)
            cur_late_sum += add_late
            cur_cumul = cumul_new
            prev = idx
            cur_duration_if_close = dur_if_close_A
            # 同步簿记
            cur_payload_kg = new_payload
            cur_path_km_to_last = new_path_km_to_last
            cur_energy_payload_num = new_energy_payload_num
        else:
            # 先收掉当前趟
            # 计算当前趟的电量
            route_km = cur_path_km_to_last + tt.c2dep_dist[prev]
            energy_num = ep.curb_weight_kg * route_km + cur_energy_payload_num
            energy_frac = energy_num / denom_energy
            
            total_late_cost += ep.late_weight * cur_late_sum
            total_duration  += cur_duration_if_close
            total_energy += energy_frac
            routes.append(cur_route)
            route_energies.append(energy_frac)

            # 开新趟跑该单
            cur_route = [cid]
            cur_payload_kg = demand[idx]
            cur_path_km_to_last = tt.dep2c_dist[idx]
            cur_energy_payload_num = demand[idx] * cur_path_km_to_last

            cur_cumul = tt.dep2c[idx]
            cur_late_sum = max(0.0, cur_cumul - tt.deadline[idx]) * tt.penalty_w[idx]
            cur_cumul += tt.service[idx]
            prev = idx
            
            # 计算新趟电量和充电时间
            route_km = cur_path_km_to_last + tt.c2dep_dist[prev]
            energy_num = ep.curb_weight_kg * route_km + cur_energy_payload_num
            energy_frac = energy_num / denom_energy
            charge_time = energy_frac / 0.04
            cur_duration_if_close = cur_cumul + tt.c2dep[prev] + charge_time

            if (cur_payload_kg > ep.max_payload_kg) or (energy_frac > 1.0):
                return DecodeResult(INF, INF, 1e9, INF, [])

    # 收尾
    if cur_route:
        # 计算最后一趟的电量
        route_km = cur_path_km_to_last + tt.c2dep_dist[prev]
        energy_num = ep.curb_weight_kg * route_km + cur_energy_payload_num
        energy_frac = energy_num / denom_energy
        
        total_late_cost += ep.late_weight * cur_late_sum
        total_duration  += cur_duration_if_close
        total_energy += energy_frac
        routes.append(cur_route)
        route_energies.append(energy_frac)

    overtime = max(0.0, total_duration - work_budget)
    total_cost = total_late_cost + ep.lambda_duration * total_duration + ep.overtime_penalty * overtime
    return DecodeResult(total_cost, total_late_cost, total_duration, total_energy, routes, route_energies)

# =======================
# 评估：90% 快评估 + 10% 精评（每代最优）
# =======================
USE_GREEDY_RATIO = 0.9

def eval_solution(sol: Solution, id2idx_map, tt: TravelTime, ep: EvalParams,
                  work_budget: float, force_dp=False) -> DecodeResult:
    if (not force_dp) and (random.random() < USE_GREEDY_RATIO):
        return split_greedy(sol.order, id2idx_map, tt, ep, work_budget)
    else:
        return split_dp(sol.order, id2idx_map, tt, ep, work_budget)

# =======================
# GA + Tabu（隔代 Tabu & First-Improvement）
# =======================
def init_solution(n):
    arr = list(range(1, n+1))
    random.shuffle(arr)
    return Solution(arr)

def order_crossover(p1: Solution, p2: Solution) -> Solution:
    n = len(p1.order)
    a, b = sorted(random.sample(range(n), 2))
    child = [None]*n
    child[a:b+1] = p1.order[a:b+1]
    fill = [g for g in p2.order if g not in child]
    k = 0
    for i in range(n):
        if child[i] is None:
            child[i] = fill[k]; k += 1
    return Solution(child)

def mutate(sol: Solution, swap_rate=0.2, shuffle_block_rate=0.05):
    arr = sol.order
    nswap = max(1, int(len(arr) * swap_rate * random.random()))
    for _ in range(nswap):
        i, j = random.sample(range(len(arr)), 2)
        arr[i], arr[j] = arr[j], arr[i]
    if random.random() < shuffle_block_rate and len(arr) >= 6:
        a, b = sorted(random.sample(range(len(arr)), 2))
        sub = arr[a:b+1]
        random.shuffle(sub)
        arr[a:b+1] = sub

def tournament_select(popz, k=3) -> Solution:
    cand = random.sample(popz, k)
    cand.sort(key=lambda x: x[1].cost)
    return copy.deepcopy(cand[0][0])

def tabu_improve(sol: Solution, id2idx_map, tt: TravelTime, ep: EvalParams, work_budget: float,
                 tenure=8, iters=25, neigh_k=20, neigh_window: int = None):
    best_sol = copy.deepcopy(sol)
    best_eval = eval_solution(best_sol, id2idx_map, tt, ep, work_budget, force_dp=False)
    cur_sol, cur_eval = copy.deepcopy(best_sol), copy.deepcopy(best_eval)
    tabu = {}  # move -> expire_step

    def add_tabu(mv, step):
        tabu[mv] = step + tenure

    step = 0
    n = len(cur_sol.order)

    while step < iters:
        step += 1
        moved = False

        tries = 0
        while tries < neigh_k:
            tries += 1
            if random.random() < 0.5:
                # swap
                if neigh_window:
                    i = random.randrange(n)
                    j = min(n-1, max(0, i + random.randint(-neigh_window, neigh_window)))
                    if i == j: continue
                    i, j = sorted((i, j))
                else:
                    i, j = sorted(random.sample(range(n), 2))
                mv = ("s", cur_sol.order[i], cur_sol.order[j])
                ns = copy.deepcopy(cur_sol)
                ns.order[i], ns.order[j] = ns.order[j], ns.order[i]
            else:
                # relocate
                i = random.randrange(n)
                if neigh_window:
                    j = min(n-1, max(0, i + random.randint(-neigh_window, neigh_window)))
                    if i == j: continue
                else:
                    j = random.randrange(n)
                    if i == j: continue
                mv = ("r", cur_sol.order[i], i, j)
                ns = copy.deepcopy(cur_sol)
                gene = ns.order.pop(i)
                ns.order.insert(j, gene)

            ev = eval_solution(ns, id2idx_map, tt, ep, work_budget, force_dp=False)
            if (tabu.get(mv, -1) >= step) and (ev.cost >= best_eval.cost - 1e-9):
                continue

            cur_sol, cur_eval = ns, ev
            add_tabu(mv, step)
            moved = True
            if cur_eval.cost < best_eval.cost - 1e-9:
                best_sol, best_eval = copy.deepcopy(cur_sol), copy.deepcopy(cur_eval)
            break

        if not moved:
            mutate(cur_sol, swap_rate=0.05)
            cur_eval = eval_solution(cur_sol, id2idx_map, tt, ep, work_budget, force_dp=False)

    return best_sol, best_eval

@dataclass
class GAParams:
    pop_size: int = 60
    generations: int = 2000
    crossover_rate: float = 0.9
    mutation_rate: float = 0.30
    elite_frac: float = 0.30
    tabu_elite_k: int = 1
    tabu_iters: int = 12
    tabu_tenure: int = 8
    neigh_k: int = 12
    seed: int = 123
    run_tabu_every: int = 12
    neigh_window: int = 6
    
    early_patience: int = 200          # 连续无改进代数
    early_min_delta: float = 1e-4      # 认为“有改进”的最小下降幅度
    early_start_gen: int = 50          # 多少代之后才开始计无改进（防早期误停）

# =======================
# 带退火的 solve（预算放宽→收紧；罚金升高）
# =======================
def solve(ins: Instance, ep: EvalParams, gp: GAParams):
    seed_all(gp.seed)
    tt = TravelTime(ins)
    id2idx = {c.id: i for i, c in enumerate(ins.customers)}

    # ———— 退火参数（可按需调）————
    relax_start = 1.8
    tighten_until = 0.7
    pen_start = 1e6
    pen_end   = 1e9

    def schedule(gen, N):
        phase = gen / max(1, N-1)
        if phase < tighten_until:
            factor = relax_start - (relax_start - 1.0) * (phase / tighten_until)
        else:
            factor = 1.0
        cur_budget = ins.work_budget * factor
        cur_penalty = pen_start * (pen_end / pen_start) ** phase
        return cur_budget, cur_penalty

    NGEN = gp.generations
    cur_budget, cur_penalty = schedule(0, NGEN)
    ep_cur = replace(ep, overtime_penalty=cur_penalty)

    pop = [init_solution(len(ins.customers)) for _ in range(gp.pop_size)]
    evals = [eval_solution(s, id2idx, tt, ep_cur, cur_budget, force_dp=False) for s in pop]
    popz = list(zip(pop, evals))

    bsol, beval = min(popz, key=lambda x: x[1].cost)
    beval = eval_solution(bsol, id2idx, tt, ep_cur, cur_budget, force_dp=True)
    best_sol, best_eval = copy.deepcopy(bsol), copy.deepcopy(beval)

    for gen in range(NGEN):
        cur_budget, cur_penalty = schedule(gen, NGEN)
        ep_cur = replace(ep, overtime_penalty=cur_penalty)

        popz.sort(key=lambda x: x[1].cost)
        elite_k = max(1, int(gp.pop_size * gp.elite_frac))
        new_popz = popz[:elite_k]

        if (gen % 50) == 0:
            print(f"[sched] gen={gen} budget={cur_budget:.1f} penalty={cur_penalty:.2e}")

        if (gen % gp.run_tabu_every) == 0:
            for i in range(min(gp.tabu_elite_k, elite_k)):
                e_sol = copy.deepcopy(new_popz[i][0])
                imp_sol, imp_eval = tabu_improve(
                    e_sol, id2idx, tt, ep_cur, cur_budget,
                    tenure=gp.tabu_tenure, iters=gp.tabu_iters,
                    neigh_k=gp.neigh_k, neigh_window=gp.neigh_window
                )
                imp_eval = eval_solution(imp_sol, id2idx, tt, ep_cur, cur_budget, force_dp=True)
                if imp_eval.cost < new_popz[i][1].cost - 1e-9:
                    new_popz[i] = (imp_sol, imp_eval)

        while len(new_popz) < gp.pop_size:
            p1 = tournament_select(popz, 3); p2 = tournament_select(popz, 3)
            child = order_crossover(p1, p2) if random.random() < gp.crossover_rate else copy.deepcopy(p1)
            if random.random() < gp.mutation_rate:
                mutate(child, swap_rate=0.2, shuffle_block_rate=0.05)
            ch_eval = eval_solution(child, id2idx, tt, ep_cur, cur_budget, force_dp=False)
            new_popz.append((child, ch_eval))

        popz = new_popz

        cur_best_sol, _ = min(popz, key=lambda x: x[1].cost)
        cur_best_eval = eval_solution(cur_best_sol, id2idx, tt, ep_cur, cur_budget, force_dp=True)
        if cur_best_eval.cost < best_eval.cost - 1e-9:
            best_sol, best_eval = copy.deepcopy(cur_best_sol), copy.deepcopy(cur_best_eval)

        if (gen + 1) % 10 == 0:
            print(f"[Gen {gen+1}] best={best_eval.cost:.3f} late={best_eval.total_late_cost:.3f} "
                  f"dur={best_eval.total_duration:.1f}m trips={len(best_eval.routes)}")

    return best_sol, best_eval

# =======================
# 数据构建（与你的数据清洗保持一致）
# =======================
def build_instance_from_df(df: pd.DataFrame,
                           XCOL: str, YCOL: str, LCOL: str, PCOL: str, WCOL: str,
                           depot_xy=(0.0, 0.0),
                           speed_kmh: float = 40.0,
                           work_budget_min: float = 3200.0) -> Instance:
    customers = []
    for i, row in df.iterrows():
        customers.append(Customer(
            id=i+1,
            x=float(row[XCOL]),
            y=float(row[YCOL]),
            service=0.0,
            deadline=float(row[LCOL]),
            weight=float(row[PCOL]),     # 元/分钟
            demand=float(row[WCOL])     # kg
        ))
    dep = Depot(float(depot_xy[0]), float(depot_xy[1]))
    ins = Instance(
        customers=customers,
        depot=dep,
        speed=speed_kmh/60.0,            # km/h → km/min
        work_budget=work_budget_min
    )
    return ins

def print_solution(res: DecodeResult):
    print("\n=== Best Plan (Single vehicle, multi-trips) ===")
    total_charging_time = res.total_energy / 0.04  # 总充电时间
    print(f"Total Cost: {res.cost:.3f} | LateCost: {res.total_late_cost:.3f} | "ION_ID = 0):
          f"Total Work: {res.total_duration:.1f} min | Total Energy: {res.total_energy:.2f} | "
          f"Charging Time: {total_charging_time:.1f} min | #Trips: {len(res.routes)}")
    
    for k, r in enumerate(res.routes, 1):
        path = "0 -> " + " -> ".join(map(str, r)) + " -> 0"
        energy_pct = res.route_energies[k-1] * 100 if res.route_energies else 0
        charge_time = res.route_energies[k-1] / 0.04 if res.route_energies else 0
        print(f"  Trip {k}: {path} | Energy: {energy_pct:.1f}% | Charge Time: {charge_time:.1f} min")

def save_results_to_file(res: DecodeResult, sheet_name: str, region_id: int, filename: str = None):
    """保存结果到txt文件"""
    if filename is None:
        filename = f"/Users/linshangjin/25CCM/NKU/results_{sheet_name}_region_{region_id}.txt"
    
    with open(filename, 'w', encoding='utf-8') as f:
        f.write("=== Vehicle Routing Results ===\n")
        f.write(f"Sheet: {sheet_name}, Region: {region_id}\n\n")
        
        total_charging_time = res.total_energy / 0.04
        f.write(f"Total Cost: {res.cost:.3f}\n")
        f.write(f"Late Cost: {res.total_late_cost:.3f}\n")
        f.write(f"Total Work Time: {res.total_duration:.1f} min\n")
        f.write(f"Total Energy Consumed: {res.total_energy:.2f} (equivalent to {res.total_energy*100:.1f}%)\n")
        f.write(f"Total Charging Time: {total_charging_time:.1f} min\n")
        f.write(f"Number of Trips: {len(res.routes)}\n\n")
        
        f.write("=== Trip Details ===\n")
        for k, r in enumerate(res.routes, 1):
            path = "0 -> " + " -> ".join(map(str, r)) + " -> 0"*480)
            energy_pct = res.route_energies[k-1] * 100 if res.route_energies else 0
            charge_time = res.route_energies[k-1] / 0.04 if res.route_energies else 0
            f.write(f"Trip {k}: {path}\n")
            f.write(f"  Energy Used: {energy_pct:.1f}%\n")
            f.write(f"  Charging Time: {charge_time:.1f} min\n\n")
        
        f.write("=== Summary ===\n")
        f.write(f"Total Late Cost: {res.total_late_cost:.3f}\n")
        f.write(f"Total Duration: {res.total_duration:.1f} min\n")
        f.write(f"Total Energy: {res.total_energy:.3f}\n")

def GA_TS(SHEET_NAME = "Case1", REGION_ID = 0):
    PATH_XLSX = "附件一_A_data.xlsx"
    SHEET_NAME = "Case1"
    PATH_JSON = "region_binaries.json"
    REGION_ID = 0

    with open(PATH_JSON, "r", encoding="utf-8") as f:
        region_bins = json.load(f)
    idx_mask = region_bins[f"region_{REGION_ID}"]

    df = pd.read_excel(PATH_XLSX, sheet_name=SHEET_NAME).copy()
    df.columns = df.columns.str.strip().str.replace(r"[^\w]", "_", regex=True)

    mask = [bool(x) for x in idx_mask]
    df = df.loc[mask].reset_index(drop=True)

    XCOL = "x坐标__km_"
    YCOL = "y坐标__km_"
    LCOL = "规定配送时间__分钟_"
    PCOL = "延迟处罚成本__元_分钟_"
    WCOL = "载重__kg_"

    for col in [XCOL, YCOL, LCOL, PCOL, WCOL]:
        df[col] = pd.to_numeric(df[col], errors="coerce")
    df = df.dropna(subset=[XCOL, YCOL, LCOL, PCOL, WCOL]).reset_index(drop=True)

    ins = build_instance_from_df(df, XCOL, YCOL, LCOL, PCOL, WCOL,
                                 depot_xy=(6.75, 6.75),
                                 speed_kmh=40.0,
                                 work_budget_min=3200.0)

    ep  = EvalParams(
        late_weight=1.0,
        lambda_duration=0.02,
        overtime_penalty=1e9,
        max_total_load_kg=500.0,
        curb_weight_kg=100.0,
        max_payload_kg=400.0,
        energy_max_range_km=100.0
    )

    gp = GAParams(
        pop_size=60, generations=2000,
        crossover_rate=0.9, mutation_rate=0.30,
        elite_frac=0.30,
        tabu_elite_k=1, tabu_iters=12, tabu_tenure=8,
        neigh_k=12, run_tabu_every=12, neigh_window=6,
        seed=123
    )

    best_sol, best_eval = solve(ins, ep, gp)
    print_solution(best_eval)
    
    # 保存结果到文件
    save_results_to_file(best_eval, SHEET_NAME, REGION_ID)
    print(f"\nResults saved to: results_{SHEET_NAME}_region_{REGION_ID}.txt")
    
    # 返回总迟到成本、总工作时长和总电量
    return best_eval.total_late_cost, best_eval.total_duration, best_eval.total_energy
